Shifts in bird and butterfly communities of Tamhini Ghat
Introduction
This file was last update on 2021-06-30
#Required libraries
library(tidyverse)
library(broom)
library(lubridate)
library(viridis)
library(RColorBrewer)
#Data import
abn <- read.csv("./Data/Abundance_280820.csv")
transects <- read.csv("./Data/Transects.csv")
# adding relevant variables
abn <- abn%>%
mutate(Year = year(dmy(Date)))%>% #adding year of sampling
mutate(Study = case_when(Year %in% c(1998, 1999, 2000, 2001) ~ "1998 - 2000", #previous study
Year %in% c(2016, 2017) ~ "2016 - 2017"))%>% #current study
mutate(Transect = as.character(Transect),
Transect = ifelse(Transect == "T6A" | Transect == "T6B", "T6", Transect),
Transect = as.factor(Transect))
transects <- transects%>%
mutate(Transect = as.factor(Transect))
# cleaning abundance
abn <- abn%>%
filter(!is.na(Abundance) &
Specific.Epithet != "" & #removing samples with no species
Specific.Epithet != " " &
Specific.Epithet != "Unidentified" & #removing unidentified species
substr(Specific.Epithet, 1, 1) != "?")%>% #83 obs removed
droplevels()
theme_set(theme_bw()+
theme(legend.position = "top",
text = element_text(size = 20)
)
)#setting ggplot theme
# function for t-test reporting
tstats <- function(x){
x%>%
rename('1998-2001' = estimate1,
'2016-2017' = estimate2)%>%
dplyr::select(`1998-2001`, `2016-2017`, parameter, statistic, p.value)
}Sampling effort
No. of visits to each transect/site.
abn%>%
group_by(Study, Taxa, Transect)%>%
summarise(n = length(unique(Date)))%>%
spread(Transect, n)| Study | Taxa | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| 1998 - 2000 | Birds | 19 | 18 | 7 | 10 | 4 | 1 | 1 |
| 1998 - 2000 | Butterflies | 14 | 7 | 4 | 9 | 4 | 1 | 1 |
| 2016 - 2017 | Birds | 23 | 23 | 22 | 19 | 16 | 22 | 17 |
| 2016 - 2017 | Butterflies | 23 | 20 | 20 | 22 | 18 | 21 | 15 |
No. of individuals sampled
abn%>%
group_by(Study, Taxa, Transect)%>%
summarise(n = sum(Abundance))%>%
spread(Transect, n)| Study | Taxa | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| 1998 - 2000 | Birds | 255 | 165 | 103 | 213 | 236 | 7 | 28 |
| 1998 - 2000 | Butterflies | 93 | 59 | 12 | 180 | 157 | 10 | 4 |
| 2016 - 2017 | Birds | 480 | 393 | 358 | 289 | 165 | 203 | 133 |
| 2016 - 2017 | Butterflies | 432 | 188 | 207 | 652 | 187 | 224 | 124 |
Uneveness in sampling effort across sites in studies may bias estimates of change.
Rarefying across sites
library(vegan)
# input for vegan
veg_birds <- abn%>%
filter(Taxa == "Birds")%>%
group_by(Study, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = paste(Study, "_", Transect))%>%
group_by(sample)%>%
#mutate(N = sum(n, na.rm = T))%>%
#filter(N > 0)%>%#removing samples where no species were detected
dplyr::select(sample, Specific.Epithet, n)%>%
ungroup()%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
veg_butterfly <- abn%>%
filter(Taxa == "Butterflies")%>%
group_by(Study, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = paste(Study, "_", Transect))%>%
group_by(sample)%>%
#mutate(N = sum(n, na.rm = T))%>%
#filter(N > 0)%>%#removing samples where no species were detected
dplyr::select(sample, Specific.Epithet, n)%>%
ungroup()%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
# sample info
samp <- abn%>%
left_join(transects, by = c("Transect"))%>%
dplyr::select(Study, Transect, Habitat)%>%
distinct()%>%
mutate(sample = paste(Study, "_", Transect))
## Species accumulation curves
rarebutts <- rarecurve(veg_butterfly) ## butterflies# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html
names(rarebutts) <- samp$sample
protox <- mapply(
FUN = function(x, y) {
mydf <- as.data.frame(x)
colnames(mydf) <- "value"
mydf$sites <- y
mydf$subsample <- attr(x, "Subsample")
mydf},
x = rarebutts, y = as.list(names(rarebutts)), SIMPLIFY = FALSE)
rarebutts.df <- do.call(rbind, protox)
rownames(rarebutts.df) <- NULL # pretty
rarebutts.df$Taxa <- "Butterflies"
## repeating for birds
rarebirds <- rarecurve(veg_birds)# Coerce data into "long" form.
## taken from https://stat.ethz.ch/pipermail/r-sig-ecology/2018-December/005867.html
names(rarebirds) <- samp$sample
protox <- mapply(
FUN = function(x, y) {
mydf <- as.data.frame(x)
colnames(mydf) <- "value"
mydf$sites <- y
mydf$subsample <- attr(x, "Subsample")
mydf},
x = rarebirds, y = as.list(names(rarebirds)), SIMPLIFY = FALSE)
rarebirds.df <- do.call(rbind, protox)
rownames(rarebirds.df) <- NULL # pretty
rarebirds.df$Taxa <- "Birds"
## joining rarefaction values
rarefied <- full_join(rarebirds.df, rarebutts.df)
## adding sample info
rarefied <- left_join(rarefied, samp, by = c("sites" = "sample"))# Plot.
ggplot(rarefied, aes(x = subsample, y = value, col = Study))+
geom_line(size = 1)+
labs(x = "Sample Size", y = "Species")+
facet_grid(Taxa ~ Transect)ggsave(last_plot(), filename = "./Figures and Tables/figureA.png", height = 8, width = 16, dpi = 300) Change in diversity across years
Species richness
Across the two studies
abn%>%
group_by(Taxa, Study)%>%
summarise(richness = length(unique(Specific.Epithet)),
N = sum(Abundance, na.rm = T))| Taxa | Study | richness | N |
|---|---|---|---|
| Birds | 1998 - 2000 | 70 | 1007 |
| Birds | 2016 - 2017 | 105 | 2021 |
| Butterflies | 1998 - 2000 | 45 | 515 |
| Butterflies | 2016 - 2017 | 66 | 2014 |
More species of both taxa were encountered in current study compared to the previous study. Difference in sampling effort is apparent accross two studies.
Across sites in the two studies
No. of species encountered across sites sampled in the two studies.
abn%>%
left_join(transects, by = c("Transect"))%>%
group_by(Taxa, Study, Transect)%>%
summarise(richness = length(unique(Specific.Epithet)))%>%
spread(Transect, richness)| Taxa | Study | T1A | T1B | T2 | T3 | T4 | T5 | T6 |
|---|---|---|---|---|---|---|---|---|
| Birds | 1998 - 2000 | 47 | 37 | 28 | 36 | 29 | 4 | 11 |
| Birds | 2016 - 2017 | 66 | 59 | 53 | 47 | 31 | 38 | 28 |
| Butterflies | 1998 - 2000 | 27 | 22 | 8 | 37 | 24 | 5 | 3 |
| Butterflies | 2016 - 2017 | 41 | 32 | 36 | 44 | 31 | 39 | 27 |
Fewer species were encountered in sites with lower effort in the previous study, i.e., T5 and T6.
Shannon diversity
# Function to caluculate Shannon diversity
shanon <- function(c){
## c is a vector of species abundances in a site
## Spade cannont compute n < 20
if(sum(c) < 20){
H <- "Sample size too small"
}else{
require(SpadeR)
require(tidyr)
div <- Diversity(c, datatype = "abundance")
H <- data.frame(div$Shannon_diversity)%>%
rownames_to_column("Method")
}
return(H)
}
# Calculating diversity in each site over years
div_year <- abn%>%
group_by(Study, Taxa, Transect, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
spread(Specific.Epithet, n, fill = c(0))%>%
group_by(Study, Taxa, Transect)%>%
nest()%>%
mutate(H = map(data, ~shanon(.)))%>%
dplyr::select(-data)%>%
unnest(H)%>%
filter(grepl("MLE", Method) | is.na(Method))
# Summary
div_year%>%
group_by(Taxa, Study)%>%
skimr::skim(Estimate)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 17.31 | 6.47 |
| Birds | 2016 - 2017 | 24.88 | 5.74 |
| Butterflies | 1998 - 2000 | 15.87 | 3.42 |
| Butterflies | 2016 - 2017 | 20.40 | 3.24 |
Linear model testing change in diversity across years:
div_year%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~lm(Estimate ~ Study, data = .)),
sumry = map(test, broom::tidy),
stat = map(test, broom::glance))%>%
dplyr::select(sumry, stat)%>%
unnest()%>%
select(Taxa:r.squared)## Error in (function (classes, fdef, mtable) : unable to find an inherited method for function 'select' for signature '"grouped_df"'
Plotting diversity across studies:
div_year%>%
ggplot(aes(Taxa, Estimate, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Effective number of species")+
scale_fill_brewer(palette = "Greys")ggsave(last_plot(), filename = "./Figures and Tables/figure2.png", height = 6, width = 8, dpi = 300) The first order diversity of both taxa has increased across the two surveys. SIgnificantly for birds but not for butterflies. Likely due to low sampling effort.
Change in composition of species accross years
Turn over of each taxa across studies:
# sample info
samp_birds <- abn%>%
left_join(transects, by = c("Transect"))%>%
filter(Taxa == "Birds")%>%
dplyr::select(Study, Taxa, Transect, Habitat)%>%
distinct()%>%
mutate(sample = paste(Study, "_", Transect))
samp_butterfly <- abn%>%
left_join(transects, by = c("Transect"))%>%
filter(Taxa == "Butterflies")%>%
dplyr::select(Study, Taxa, Transect, Habitat)%>%
distinct()%>%
mutate(sample = paste(Study, "_", Transect))
# Similarity matrix by studies
sim_birds <- abn%>%
filter(Taxa == "Birds")%>%
group_by(Study, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = Study)%>%
ungroup()%>%
dplyr::select(sample, Specific.Epithet, n)%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
sim_butterfly <- abn%>%
filter(Taxa == "Butterflies")%>%
group_by(Study, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
mutate(sample = Study)%>%
ungroup()%>%
dplyr::select(sample, Specific.Epithet, n)%>%
spread(Specific.Epithet, n, fill = c(n = 0))%>%
column_to_rownames(var = "sample")
# bray - curtis
data_frame(
Taxa = c("Birds", "Butterflies"),
Dissimilarity = c(as.numeric(vegdist(sim_birds, method = "bray")),
as.numeric(vegdist(sim_butterfly, method = "bray")))
)# ordination
NMDS_birds <- metaMDS(veg_birds, distance = "bray", k = 2, model = "local")
NMDS_butterfly <- metaMDS(veg_butterfly, distance = "bray", trymax = 1000, k = 2, model = "local")
# Extracting scores and adding habitat vars
ord_birds <- data.frame(scores(NMDS_birds))%>%
rownames_to_column("sample")%>%
left_join(samp_birds, "sample")
ord_butterfly <- data.frame(scores(NMDS_butterfly))%>%
rownames_to_column("sample")%>%
left_join(samp_butterfly, "sample")
ord <- bind_rows(ord_birds, ord_butterfly)Testing change in composition of birds :
adonis(veg_birds ~ Study, data = samp_birds)##
## Call:
## adonis(formula = veg_birds ~ Study, data = samp_birds)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Study 1 0.77547 0.77547 3.9632 0.24827 0.002 **
## Residuals 12 2.34798 0.19567 0.75173
## Total 13 3.12345 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Testing change in composition of butterflies:
adonis(veg_butterfly ~ Study, data = samp_butterfly)##
## Call:
## adonis(formula = veg_butterfly ~ Study, data = samp_butterfly)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## Study 1 0.9191 0.91908 4.054 0.25253 0.001 ***
## Residuals 12 2.7205 0.22671 0.74747
## Total 13 3.6396 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
NMDS ordination plot:
# plotting
ord%>%
ggplot(aes(NMDS1, NMDS2, col = Study, fill = Study))+
stat_ellipse(level = 0.95, geom = "polygon", linetype = "dashed", size = 1, alpha = 0.3)+
geom_point(aes(shape = Habitat), size = 3)+
guides(shape = guide_legend(nrow = 2, byrow = T))+
scale_color_brewer(palette = "Set1")+
scale_fill_brewer(palette = "Set1")+
facet_wrap(~Taxa)+
theme(legend.box = "vertical")ggsave(last_plot(), filename = "./Figures and Tables/figure3.png", height = 6, width = 8, dpi = 300)Both bird and butterfly communities compositions are significantly different compared to the previous study.
Change in niche width of community
Trophic Guilds
Diversity across trophic guilds:
# compiling trophic guild data
host <- read.csv("./Data/Butterflies_host plant_280820.csv") #butterfly host plant data
diet <- read.csv("./Data/Birds_diet.csv") #bird diet data
trophic <- host%>%
full_join(diet)%>%
dplyr::select(Specific.Epithet, Guild)%>%
distinct()
# Calculating diversity of guilds
guilds <- abn%>%
left_join(trophic, by = c("Specific.Epithet"))%>%
filter(!Guild == "")%>%
group_by(Study, Taxa, Transect, Guild, Specific.Epithet)%>%
summarise(n = sum(Abundance, na.rm = T))%>%
spread(Specific.Epithet, n, fill = c(0))%>%
group_by(Study, Taxa, Transect, Guild)%>%
nest()%>%
mutate(H = map(data, ~shanon(.)))%>%
dplyr::select(-data)%>%
unnest(H)%>%
filter(grepl("MLE", Method) | is.na(Method))
# summary
guilds%>%
group_by(Taxa, Guild, Study)%>%
skimr::skim(Estimate)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Guild, Study, mean, sd)Variable type: numeric
| Taxa | Guild | Study | mean | sd |
|---|---|---|---|---|
| Birds | Carnivore | 1998 - 2000 | 3.43 | 0.31 |
| Birds | Carnivore | 2016 - 2017 | 6.26 | 0.15 |
| Birds | Frugivore | 1998 - 2000 | 4.11 | 0.92 |
| Birds | Frugivore | 2016 - 2017 | 4.21 | 0.60 |
| Birds | Grainivore | 1998 - 2000 | 3.86 | 0.75 |
| Birds | Grainivore | 2016 - 2017 | 3.84 | 0.49 |
| Birds | Insectivore | 1998 - 2000 | 8.77 | 1.95 |
| Birds | Insectivore | 2016 - 2017 | 13.58 | 4.82 |
| Birds | Omnivore | 1998 - 2000 | 4.82 | 0.92 |
| Birds | Omnivore | 2016 - 2017 | 5.88 | 1.15 |
| Butterflies | Generalist | 1998 - 2000 | 3.16 | 1.48 |
| Butterflies | Generalist | 2016 - 2017 | 5.95 | 1.87 |
| Butterflies | Grass Specialist | 1998 - 2000 | 1.80 | NA |
| Butterflies | Grass Specialist | 2016 - 2017 | 2.71 | 0.71 |
| Butterflies | Herb Specialist | 1998 - 2000 | 3.58 | NA |
| Butterflies | Herb Specialist | 2016 - 2017 | 4.09 | 0.86 |
| Butterflies | Liana Specialist | 1998 - 2000 | NaN | NA |
| Butterflies | Shrub Specialist | 1998 - 2000 | 4.36 | 0.77 |
| Butterflies | Shrub Specialist | 2016 - 2017 | 4.18 | 2.10 |
| Butterflies | Tree Specialist | 1998 - 2000 | 5.92 | 3.07 |
| Butterflies | Tree Specialist | 2016 - 2017 | 6.88 | 1.63 |
Linear model testing change in guild diversity:
guilds%>%
filter(Guild != "Liana Specialist")%>% #removing due to lack of obs
complete(nesting(Study, Taxa), Transect, Guild, fill = list(H = 0))%>%
group_by(Taxa, Guild)%>%
nest()%>%
mutate(test = map(data, ~lm(Estimate ~ Study, data = .)),
sumry = map(test, broom::tidy),
stat = map(test, broom::glance))%>%
dplyr::select(sumry, stat)%>%
unnest()%>%
dplyr::select(Taxa:r.squared)| Taxa | Guild | term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|---|---|
| Birds | Carnivore | (Intercept) | 3.4273333 | 0.0482458 | 71.0390510 | 0.0000000 | 0.9745403 |
| Birds | Carnivore | Study2016 - 2017 | 2.8301667 | 0.0849454 | 33.3174666 | 0.0000000 | 0.9745403 |
| Birds | Frugivore | (Intercept) | 4.1146000 | 0.1159192 | 35.4954216 | 0.0000000 | 0.0049246 |
| Birds | Frugivore | Study2016 - 2017 | 0.0966857 | 0.1517738 | 0.6370381 | 0.5258746 | 0.0049246 |
| Birds | Grainivore | (Intercept) | 3.8566667 | 0.1170495 | 32.9490144 | 0.0000000 | 0.0002457 |
| Birds | Grainivore | Study2016 - 2017 | -0.0166667 | 0.1850716 | -0.0900553 | 0.9287877 | 0.0002457 |
| Birds | Insectivore | (Intercept) | 8.7660000 | 0.6795468 | 12.8997745 | 0.0000000 | 0.2877674 |
| Birds | Insectivore | Study2016 - 2017 | 4.8127143 | 0.8628491 | 5.5777010 | 0.0000003 | 0.2877674 |
| Birds | Omnivore | (Intercept) | 4.8210000 | 0.1665640 | 28.9438346 | 0.0000000 | 0.2223562 |
| Birds | Omnivore | Study2016 - 2017 | 1.0560000 | 0.2180834 | 4.8421836 | 0.0000060 | 0.2223562 |
| Butterflies | Generalist | (Intercept) | 3.1625000 | 0.5157184 | 6.1322230 | 0.0000001 | 0.3200482 |
| Butterflies | Generalist | Study2016 - 2017 | 2.7838333 | 0.5738384 | 4.8512495 | 0.0000124 | 0.3200482 |
| Butterflies | Grass Specialist | (Intercept) | 1.7980000 | 0.2732411 | 6.5802684 | 0.0000001 | 0.2042940 |
| Butterflies | Grass Specialist | Study2016 - 2017 | 0.9124000 | 0.2921071 | 3.1235123 | 0.0034125 | 0.2042940 |
| Butterflies | Herb Specialist | (Intercept) | 3.5810000 | 0.3907790 | 9.1637464 | 0.0000000 | 0.0300820 |
| Butterflies | Herb Specialist | Study2016 - 2017 | 0.5111429 | 0.4064163 | 1.2576829 | 0.2142349 | 0.0300820 |
| Butterflies | Shrub Specialist | (Intercept) | 4.3575000 | 0.4412280 | 9.8758466 | 0.0000000 | 0.0024700 |
| Butterflies | Shrub Specialist | Study2016 - 2017 | -0.1781000 | 0.5220680 | -0.3411433 | 0.7345161 | 0.0024700 |
| Butterflies | Tree Specialist | (Intercept) | 5.9220000 | 0.4125034 | 14.3562448 | 0.0000000 | 0.0525717 |
| Butterflies | Tree Specialist | Study2016 - 2017 | 0.9577143 | 0.4930359 | 1.9424840 | 0.0562211 | 0.0525717 |
Plotting guild diversity across years:
guilds%>%
ggplot(aes(reorder(Guild, Estimate), Estimate, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5, position = position_dodge(preserve = "single"))+
labs(y = "Effective number of species", x = "Trophic Guild")+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
scale_fill_brewer(palette = "Greys")+
facet_wrap(~Taxa, scales = "free_x", ncol = 1)ggsave(last_plot(), filename = "./Figures and Tables/figure4.png", height = 9, width = 8, dpi = 300) Carnivorous, Grainivorous birds have reduced in the comunity, compared to insectivorous birds which have increased. Similarly, Shrub specialist and tree specialist butterflies have reduced where grass specialists and herb specialist have increased. Suprsingly, generalist species of both taxa appear to be unchanged.
Change in niche specialisation
-Can we incorporate abunndance into community specialisation index?
Trophic niche
Trophic specialisation based on Juliard et al. 2006
#
butterflies_TSI <- host%>%
mutate(H = length(unique(Hostplant.Family)))%>%
group_by(Specific.Epithet)%>%
mutate(h = length(unique(Hostplant.Family)))%>%
ungroup()%>%
mutate(TSI = sqrt((H/h)-1))
birds_TSI <- diet%>%
mutate(H = length(unique(Prey)))%>%
group_by(Specific.Epithet)%>%
mutate(h = length(unique(Prey)))%>%
ungroup()%>%
mutate(TSI = sqrt((H/h)-1))
CSI.T <- abn%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
group_by(Study, Taxa, Transect)%>%
distinct()%>%
left_join(butterflies_TSI, by = c("Specific.Epithet"))%>%
left_join(birds_TSI, by = c("Specific.Epithet"))%>%
mutate(TSI = ifelse(is.na(TSI.x), TSI.y, TSI.x))%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet, TSI)%>%
summarise(CSI.T = mean(TSI, na.rm = T))
# summary
CSI.T%>%
group_by(Taxa, Study)%>%
skimr::skim(CSI.T)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 2.59 | 0.16 |
| Birds | 2016 - 2017 | 2.53 | 0.09 |
| Butterflies | 1998 - 2000 | 5.34 | 0.39 |
| Butterflies | 2016 - 2017 | 5.07 | 0.24 |
Linear model testing difference is trophic specialization:
CSI.T%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~lm(CSI.T ~ Study, data = .)),
sumry = map(test, broom::tidy),
stat = map(test, broom::glance))%>%
dplyr::select(sumry, stat)%>%
unnest()%>%
dplyr::select(Taxa:r.squared)| Taxa | term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|---|
| Birds | (Intercept) | 2.5880088 | 0.0493181 | 52.4758739 | 0.0000000 | 0.051883 |
| Birds | Study2016 - 2017 | -0.0565188 | 0.0697463 | -0.8103492 | 0.4335164 | 0.051883 |
| Butterflies | (Intercept) | 5.3385626 | 0.1229059 | 43.4361931 | 0.0000000 | 0.170905 |
| Butterflies | Study2016 - 2017 | -0.2733716 | 0.1738151 | -1.5727722 | 0.1417534 | 0.170905 |
CSI.T%>%
ggplot(aes(Taxa, CSI.T, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Community specialisation index \n (Trophic)")+
scale_fill_brewer(palette = "Greys")+
scale_y_sqrt()Niche width of the community did not change over the two studies.
Habitat Niche
HSI <- abn%>%
left_join(transects, by = c("Transect"))%>%
group_by(Study, Taxa, Habitat, Specific.Epithet)%>%
# calculating abundance of each species in each habitat type
summarise(n = sum(Abundance, na.rm = T))%>%
group_by(Study, Taxa)%>%
mutate(N = sum(n), # total number of indviduals encountered
rel.prop = n/N)%>% # relative proportion of species
group_by(Study, Specific.Epithet)%>%
summarise(HSI = sd(rel.prop)/mean(rel.prop))%>%
replace_na(replace = list(HSI = 0))
CSI.H <- abn%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet)%>%
group_by(Study, Taxa, Transect)%>%
distinct()%>%
left_join(HSI, by = c("Study", "Specific.Epithet"))%>%
dplyr::select(Study, Transect, Taxa, Specific.Epithet, HSI)%>%
summarise(CSI.H = mean(HSI, na.rm = T))
# summary
CSI.H%>%
group_by(Taxa, Study)%>%
skimr::skim(CSI.H)%>%
skimr::yank("numeric")%>%
dplyr::select(Taxa, Study, mean, sd)Variable type: numeric
| Taxa | Study | mean | sd |
|---|---|---|---|
| Birds | 1998 - 2000 | 0.54 | 0.10 |
| Birds | 2016 - 2017 | 0.49 | 0.02 |
| Butterflies | 1998 - 2000 | 0.61 | 0.09 |
| Butterflies | 2016 - 2017 | 0.57 | 0.03 |
Linear model testing for habitat specialisation:
CSI.H%>%
group_by(Taxa)%>%
nest()%>%
mutate(test = map(data, ~lm(CSI.H ~ Study, data = .)),
sumry = map(test, broom::tidy),
stat = map(test, broom::glance))%>%
dplyr::select(sumry, stat)%>%
unnest()%>%
dplyr::select(Taxa:r.squared)| Taxa | term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|---|
| Birds | (Intercept) | 0.5363977 | 0.0278655 | 19.249520 | 0.0000000 | 0.0869142 |
| Birds | Study2016 - 2017 | -0.0421175 | 0.0394078 | -1.068760 | 0.3062059 | 0.0869142 |
| Butterflies | (Intercept) | 0.6110527 | 0.0244089 | 25.033991 | 0.0000000 | 0.1087676 |
| Butterflies | Study2016 - 2017 | -0.0417742 | 0.0345194 | -1.210166 | 0.2495141 | 0.1087676 |
CSI.H%>%
ggplot(aes(Taxa, CSI.H, fill = Study))+
geom_boxplot(width = 0.25, alpha = 0.5)+
labs(y = "Community specialisation index \n (Habitat)")+
scale_fill_brewer(palette = "Greys")Habitat niche of the community did not change over the study periods.
Seasonal changes in relative abundance
season <- abn%>%
mutate(Month = month(dmy(Date), label = T))%>%
group_by(Study, Month, Taxa)%>%
summarise(n = sum(Abundance, na.rm=T))%>%
group_by(Study, Taxa)%>%
mutate(N = sum(n),
rel.prop = n/N)
season%>%
ggplot(aes(Month, rel.prop, col = Taxa, shape = Taxa))+
geom_point(size = 3)+
geom_path(aes(group = Taxa), size = 1, linetype = "dashed")+
scale_color_brewer(palette = "Set1")+
facet_wrap(~Study, ncol = 1)Site Map
library(raster)
library(sf)
# Raw point data
points <- read.csv("./Data/sites.csv")%>%
drop_na()%>%
dplyr::select(-Location)
# Converting to SF points
points_sf <- st_as_sf(points, coords = c('Lon', 'Lat'))
st_crs(points_sf) <- 4326
# Making transect lines
lines <- points_sf%>%
mutate(Transect = as.factor(Transect))%>%
group_by(Transect)%>%
summarise()%>%
st_cast("LINESTRING")
# getting study extent
study_ext <- st_bbox(points_sf)
ext <- st_as_sfc(study_ext)
# getting osm map
library(osmdata)
# Tamhini WLS
tamhini_raw <- opq(bbox = study_ext, timeout = 100)%>%
add_osm_feature(key = 'boundary', value = 'protected_area')%>%
osmdata_sf()%>%
unique_osmdata()
tamhini <- tamhini_raw$osm_multipolygons
st_crs(tamhini) <- "+init=epsg:4326"
# Roads
tamhini_area <- st_bbox(tamhini)
roads_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'highway')%>%
osmdata_sf()%>%
unique_osmdata()
roads <- roads_raw$osm_lines
st_crs(roads) <- 4326
# Landuse
landuse_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'natural')%>%
osmdata_sf()%>%
unique_osmdata()
forest <- landuse_raw$osm_polygons%>%
filter(natural %in% c("wood", "tree"))
wood <- landuse_raw$osm_multipolygons%>%
filter(natural == "wood")
water <- landuse_raw$osm_multipolygons%>%
filter(natural == "water")
water_small <- landuse_raw$osm_polygons%>%
filter(natural == "water")
st_crs(forest) <- 4326
# waterways
waterway_raw <- opq(bbox = tamhini_area, timeout = 100)%>%
add_osm_feature(key = 'waterway')%>%
osmdata_sf()%>%
unique_osmdata()
waterways <- waterway_raw$osm_lines
st_crs(waterways) <- 4326
# elevation
library(elevatr)
library(rgdal)
tam_elev <- get_elev_raster(locations = tamhini, z = 8)
tam_contours <- rasterToContour(tam_elev, levels = seq(0, 1000, 100))# plotting site map
library(tmap)
library(cowplot) # for plot_grid() function - good to arrange multiple plots into a grid
library(grid)
# Study site
tm_shape(forest)+ #forested lands
tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(wood)+ #forested lands
tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
tm_lines(col = "white", lty = "dotted")+
tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
tm_polygons(col = "#77A1CB", border.col = NULL)+
tm_shape(waterways)+
tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
tm_lines(alpha = 0.5, lwd = 1.5)+
tm_text("ref", size = 0.3, remove.overlap = T)+
tm_shape(tamhini)+
tm_borders(lwd = 2, col = "#165916", lty = "dashed")+
tm_text("name", size = 0.5)+
tm_shape(lines, is.master = T, bbox = study_ext+c(-0.005, -0.005,0.005,0.005))+
tm_lines(col = "red", lwd = 2)+
tm_text("Transect", size = 0.5, auto.placement = T, remove.overlap = T)+
tm_graticules(n.x = 5, lines = F, labels.inside.frame = T)+
tm_scale_bar(position = c("left", "top"), text.size = 0.5, size = 0.5) +
tm_compass(position = c('right', 'top'), size = 2)+
tm_layout(scale = 1.5, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1.5)fig1a <- grid.grab()# Tamhini WLS
tm_shape(forest)+ #forested lands
tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(tam_contours)+
tm_lines(col = "white", lty = "dotted")+
tm_text(text = "level", auto.placement = T)+
tm_shape(water)+
tm_polygons(col = "#77A1CB", border.col = NULL)+
tm_text(text = "name", remove.overlap = F, size = 0.5)+
tm_shape(water_small)+
tm_polygons(col = "#77A1CB", border.col = NULL)+
tm_text("name", remove.overlap = F)+
tm_shape(waterways)+
tm_lines(col = "#77A1CB", lwd = 2)+
tm_shape(roads)+
tm_lines(alpha = 0.5, lwd = 1.5)+
tm_shape(tamhini, is.master = T, bbox = tamhini_area+c(-0.01, -0.01,0.01,0.01))+
tm_borders(col = "#165916", lty = "dashed")+
tm_text("name", size = 0.5)+
tm_shape(ext)+
tm_borders(col = "red")+
tm_scale_bar(position = c("left", "top"), text.size = 0.5) +
tm_layout(scale = 1, legend.position = c("left", "bottom"), bg.color = "#ffdb9c", asp = 1)fig1b <- grid.grab()# legend
tm_shape(forest)+
tm_polygons(col = "#81AF81", alpha = 0.5, border.col = NULL)+
tm_shape(ext)+
tm_borders(col = "red")+
tm_add_legend(type = "fill",col = c("#81AF81", "#77A1CB", "#ffdb9c"),
labels = c("Forest", "Water", "Barren Land"), title = "Legend")+
tm_add_legend(type = "line", col = c("black", "#165916", "#77A1CB", "red"),
labels = c("Roads", "WLS boundary", "Waterway", "Transect"),
lwd = c(3,3,3,3),
lty = c(1, 2, 1, 1))+
tm_layout(legend.only = T, asp = 1, legend.position = c(0.3,0.2), scale = 2)legend <- grid::grid.grab()# making inset map
tamext <- st_as_sfc(tamhini_area + c(-.5,-.5, .5, .5))
india <- read_sf("./Data/India-States.shp")
tmap_options(check.and.fix = TRUE)
tm_shape(india)+
tm_polygons(col = "#ffdb9c")+
tm_shape(tamext)+
tm_borders(lwd = 1, col = "red")+
tm_layout(scale = 1.5, bg.color = "#77A1CB", asp = 1)fig1c <- grid.grab()# making final map
library(gridExtra)
fig1 <- grid.arrange(fig1a, fig1b, fig1c, legend, layout_matrix = rbind(c(1, 1, 1),
c(1, 1, 1),
c(2, 3, 4)))ggsave(fig1, filename = "./Figures and Tables/figure1.png", height = 8, width = 8)